Interval Truth Model

Empirical Example: Verbal Quantifiers

Authors
Affiliation

Matthias Kloft

University of Marburg

Björn S. Siepe

University of Marburg

Daniel W. Heck

University of Marburg

Published

July 18, 2025

Show the code
# Libraries
packages <- c(
  "quarto",
  "tidyverse",
  "yardstick",
  "readxl",
  "cmdstanr",
  "bayesplot",
  "posterior",
  "bayestestR",
  "rmarkdown",
  "cowplot",
  "svglite",
  "psych",
  "here",
  "osfr",
  "ggpubr",
  "ggokabeito",
  "ggh4x",
  "ggtext"
)
# in case cmdstanr is not installed
# install.packages("cmdstanr",
#                  repos = c("https://mc-stan.org/r-packages/", getOption("repos")))

if (!require("pacman")) install.packages("pacman") 
pacman::p_load(packages, update = F, character.only = T)
source("00_functions.R")

1 Load Data

Show the code
# Download cleaned data (Kloft & Heck, 2024) from OSF
if (!file.exists(here("data", "df_master_long.rds"))) {
  osf_retrieve_file("https://osf.io/7azbr") |>
    osf_download(path = here("data"))
}

# Load Data and retain all rows with the verbal quantifier task
df_long <- read_rds(here("data", "df_master_long.rds")) |>
  dplyr::filter(domain == "prob_phrase") |> 
  # recompute indices
  mutate(jj = as.integer(factor(jj)),
         ii = as.integer(factor(ii)))

df_long$name_en[df_long$name_en == "Fifty-fifty chance"] <- "fifty-fifty chance"

# item names
item_names_quantifier <- df_long %>%
  dplyr::select(jj, name_en) %>%
  distinct() %>%
  arrange(jj) %>%
  pull(name_en)

1.1 Stan Data

Show the code
# rescale the responses to get rid of zeros
df_long <- df_long %>%
  dplyr::mutate(
    x_L_rescaled = (x_L / 100 + .01) / 1.03,
    x_U_rescaled = (x_U / 100 + .02) / 1.03
  )

### Stan data declaration
df_long <- df_long %>% dplyr::filter(!is.na(x_splx_1))

I = length(unique(df_long$ii))
J = length(unique(df_long$jj))
N  = nrow(df_long)
ii = df_long$ii
jj = df_long$jj
nn = c(1:N)
Y_splx = cbind(
  df_long$x_L_rescaled,
  df_long$x_U_rescaled - df_long$x_L_rescaled,
  1 - df_long$x_U_rescaled) |> 
  as.matrix()

#Y_splx[is.na(Y_splx)]

# # sanity check: all rows sum to 1
# check <- apply(Y_splx, 1, sum) %>% as.data.frame()
# table(check)

## stan data list
stan_data <- list(
  I = I,
  J = J,
  N = N,
  ii = ii,
  jj = jj,
  nn = nn,
  Y_splx = Y_splx,
  padding = .005
)

# Choose model to fit
model_name <- "itm_beta"

2 Fit Full Model

We first fit the full model to the verbal quantifier data to see what the parameter estimates look like. We can subsequently customize the model.

2.1 Compile Model

Show the code
# Compile model
model <-
  cmdstanr::cmdstan_model(
    stan_file = here("src", "models", paste0(model_name, ".stan")),
    pedantic = TRUE,
    quiet = FALSE
  )

2.2 Fit Model

Show the code
# number of MCMC chains
n_chains <- 4
# Run sampler
fit <- model$sample(
  data = stan_data,
  seed = 2023,
  chains = n_chains,
  parallel_chains = n_chains,
  iter_warmup = 500,
  iter_sampling = 500,
  refresh = 500,
  thin = 1,
  adapt_delta = .8,
  init = .1
)

# save fit
fit$save_object(file =  here("fits", paste0(model_name, "quantifier_full_fit.RDS")))
Show the code
# load  fit
fit <- readRDS(file =  here("fits", paste0(model_name, "quantifier_full_fit.RDS")))

Show the code
parameters <- c(
  "Tr_loc",
  "Tr_wid",
  "E_loc",
  "E_wid",
  "a_loc",
  "b_loc",
  "lambda_loc",
  "lambda_wid",
  "omega",
  "mu_E",
  "sigma_I",
  "sigma_lambda",
  "rho_lambda",
  "rho_E",
  "rho_lambda"
)
Show the code
estimates_summary <- fit$summary(variables = parameters)

2.2.1 Sampler Diagnostics

Show the code
# sampler diagnostics
fit$sampler_diagnostics(format = "df") %>% 
  psych::describe(quant = c(.05,.95),) %>%
  round(2) %>%  
  as.data.frame() %>% 
  dplyr::select(median, min, Q0.05, Q0.95,  max) %>% 
  .[-c(7:9),]
              median    min  Q0.05  Q0.95     max
treedepth__     6.00   6.00   6.00   7.00    7.00
divergent__     0.00   0.00   0.00   0.00    0.00
energy__      844.95 700.31 779.25 915.71 1000.25
accept_stat__   0.94   0.41   0.70   1.00    1.00
stepsize__      0.06   0.05   0.05   0.06    0.06
n_leapfrog__   63.00  63.00  63.00 127.00  255.00
Show the code
# convergence diagnostics
convergence_summary <- 
  fit$draws(format = "df", variables = parameters) %>%
  summarise_draws(.x = ., "rhat", "ess_bulk", "ess_tail") %>%
  remove_missing() %>%
  select(-variable) %>%
  psych::describe(., quant = c(.05, .95)) %>%
  as.data.frame() %>%
  select(median, Q0.05, Q0.95, min, max)

convergence_summary %>% round(3)
           median   Q0.05    Q0.95     min      max
rhat        1.003   0.999    1.012   0.999    1.026
ess_bulk  863.624 259.573 2906.932 160.787 3952.455
ess_tail 1277.749 523.688 1641.116 344.311 1825.163

2.2.2 Effective sample size (ESS) & Rhat Plots

Show the code
# color scheme
color_scheme_set(scheme = "purple")

# Effective sample sizes
plot_neff <-
  mcmc_neff_hist(bayesplot::neff_ratio(fit, pars = parameters), binwidth = .01) +
  labs(title = "A") +
  guides(color = FALSE, fill = FALSE) +
  theme(
    legend.text = element_blank(),
    legend.key = element_blank(),
    title = element_text(size = 16, face = "bold")
  )
# Rhat
plot_rhat <-
  bayesplot::mcmc_rhat_hist(bayesplot::rhat(fit, pars = parameters)) +
  labs(title = "B") +
  guides(color = FALSE, fill = FALSE) +
  theme(
    legend.text = element_blank(),
    legend.key = element_blank(),
    title = element_text(size = 16, face = "bold")
  ) +
  yaxis_text(on = TRUE)

# Combined plot
plot_diagnostics <- gridExtra::grid.arrange(plot_neff, plot_rhat, ncol = 2)

2.3 Parameter Plots

Show the code
# color scheme
color_scheme_set(scheme = "purple")

2.3.1 Person Competence: Location

Show the code
# order estimates by size
E_loc_ordered <-
  fit$summary("E_loc") %>% 
  as.data.frame() %>% 
  arrange(median) %>% 
  select(variable) %>% 
  unlist()

# plot
mcmc_intervals(fit$draws(E_loc_ordered), point_est = "median", transformations = "log") +
  labs(
    subtitle = "Person Competence: Location",
    x = expression(E[loc]),
    y = "Respondent"
  ) +
  scale_y_discrete(labels = NULL, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

2.3.2 Person Competence: Width

Show the code
# order estimates by size of theta
E_wid_ordered <- 
  str_replace(E_loc_ordered, "E_loc", "E_wid")

# plot
mcmc_intervals(fit$draws(E_wid_ordered), point_est = "median",transformations = "log") +
  labs(subtitle = "Person Competence: Width",
       x = expression(E[wid]),
       y = "Respondent") +
  scale_y_discrete(labels = NULL, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

2.3.3 Correlation: Person Competence Location vs. Width

Show the code
mcmc_intervals(fit$draws("rho_E"), point_est = "median") +
  bayesplot::vline_0(linetype = "dashed", color = "grey70") +
  labs(
    subtitle = "Correlation: Person Competence Location vs. Width",
    x = expression(rho[E]),
    y = "Respondent"
  ) +
  scale_x_continuous(limits = c(-1,1), expand = expansion()) +
  theme_icm(base_size = 16, hide_axis_text_y = TRUE) +
  theme(panel.grid = element_blank())

2.3.4 Person Scaling Bias: Location

Show the code
# order estimates by size
a_loc_ordered <-
  fit$summary("a_loc") %>% 
  as.data.frame() %>% 
  arrange(median) %>% 
  select(variable) %>% 
  unlist()

# plot
  mcmc_intervals(fit$draws(a_loc_ordered), point_est = "median", transformations = "log") +
  labs(
    subtitle = "Person Scaling Bias: Location",
    x = expression(a[i]),
    y = "Respondent"
  ) +
  scale_y_discrete(labels = NULL, expand = expansion(.02)) +
  theme_icm(base_size = 16, hide_axis_text_y = TRUE) +
  theme(panel.grid = element_blank())

2.3.5 Person Shifting Bias: Location

Show the code
# order estimates by size
b_loc_ordered <-
  fit$summary("b_loc") %>% 
  as.data.frame() %>% arrange(median) %>% 
  select(variable) %>% 
  unlist()

# plot
  mcmc_intervals(fit$draws(b_loc_ordered), point_est = "median") +
  labs(
    subtitle = "Person Shifting Bias: Location",
    x = expression(b_loc[i]),
    y = "Respondent"
  ) +
  scale_y_discrete(labels = NULL, expand = expansion(.02)) +
  theme_icm(base_size = 16, hide_axis_text_y = TRUE) +
  theme(panel.grid = element_blank())

We find no relevant variance for the person shifting bias. The model should therefore be simplfied by ommitting the person shifting bias.

2.3.6 Person Shifting Bias: Width

Show the code
# order estimates by size of theta
b_wid_ordered <- 
  str_replace(b_loc_ordered, "b_loc", "b_wid")

# plot
mcmc_intervals(fit$draws(b_wid_ordered), point_est = "median") +
  labs(
    subtitle = "Person Shifting Bias: Width",
    x = expression(b_wid[i]),
    y = "Respondent"
  ) +
  scale_y_discrete(labels = NULL, expand = expansion(.02)) +
  theme_icm(base_size = 16, hide_axis_text_y = TRUE) +
  theme(panel.grid = element_blank())

2.3.7 Item Discernability: Location

Show the code
mcmc_intervals(fit$draws("lambda_loc"), point_est = "median", transformations = "log") +
  labs(
    subtitle = "Item Discernability: Location",
    x = expression(lambda[loc]),
    y = "Respondent"
  ) +
  scale_y_discrete(labels = item_names_quantifier, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

2.3.8 Item Discernability: Width

Show the code
mcmc_intervals(fit$draws("lambda_wid"), point_est = "median",transformations = "log") +
  labs(
    subtitle = "Item Discernability: Width",
    x = expression(lambda[wid]),
    y = "Respondent"
  ) +
  scale_y_discrete(labels = item_names_quantifier, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

2.3.9 Correlation: Iem Discernibility, Location vs. Width

Show the code
mcmc_intervals(fit$draws("rho_lambda"), point_est = "median") +
  bayesplot::vline_0(linetype = "dashed", color = "grey70") +
  labs(x = expression(rho[lambda]), y = NULL) +
  scale_x_continuous(limits = c(-1, 1), expand = expansion()) +
  theme_icm(base_size = 16, hide_axis_text_y = TRUE) +
  theme(panel.grid = element_blank())

2.3.10 Correlation of Residuals: Location vs. Width

Show the code
mcmc_intervals(fit$draws(c("omega")), point_est = "median") +
  
  labs(
    title = "B",
    subtitle = "DDRM: Expansion",
    x = expression(omega),
    y = "Respondent"
  ) +
  scale_x_continuous(limits = c(-1,1), expand = expansion()) +
  scale_y_discrete(labels = item_names_quantifier, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

2.3.11 Latent Consensus Intervals

Show the code
consensus <- data.frame(
  idx = 1:J,
  name = item_names_quantifier,
  lower = fit$summary("Tr_splx") %>% as.data.frame() %>% pull(median) %>% .[1:J],
  wid = fit$summary("Tr_splx") %>% as.data.frame() %>% pull(median) %>% .[(J +
                                                                             1):(J * 2)]
) %>%
  mutate(loc = lower + wid / 2, upper = lower + wid) %>%
  arrange(lower)

consensus %>% 
  ggplot(aes(x = loc, y = 1:J)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.3) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = seq(0, 1, .25),
    expand = expansion()
  ) +
  scale_y_continuous(breaks = 1:J, labels = item_names_quantifier[consensus$idx]) +
  labs(x = "Latent Truth",
       y = "Item") +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank(),
        axis.line.y = element_blank(),
    axis.ticks.x = element_line(colour = "#6d6d6e")
    )

Show the code
# save plot
ggsave(
  here("plots", "empirical_example", "consensus_intervals_full_model.pdf"),
  width = 15,
  height = 8,
  units = "cm",
  scale = 1.3
)
Show the code
consensus %>% 
  mutate(across(c(lower, wid, loc, upper), ~round(., 3)*100)) %>% 
  select(-idx,-wid)
                 name lower  loc upper
1               never   0.6  1.4   2.2
2        almost never   2.7  7.4  12.1
3              hardly   4.7 11.5  18.2
4              rarely   5.0 12.6  20.1
5        now and then  18.4 29.0  39.5
6        occasionally  18.7 29.6  40.6
7            possibly  19.7 33.9  48.1
8           sometimes  22.7 34.6  46.6
9         potentially  28.6 46.2  63.9
10 fifty-fifty chance  48.4 50.1  51.7
11         frequently  63.8 75.9  88.0
12              often  66.5 77.6  88.7
13   most of the time  70.3 81.2  92.1
14         constantly  77.9 87.1  96.4
15      almost always  84.4 90.8  97.2
16             always  97.0 98.2  99.4

3 Fit Customized Model

We removed the person shifting bias from the model.

3.0.1 Compile Model

Specify model name:

Show the code
# Choose model to fit
model_name_quantifier <- "itm_quantifier_beta"
Show the code
# Compile model
model_quantifier <-
  cmdstanr::cmdstan_model(
    stan_file = here("src", "models", paste0(model_name_quantifier, ".stan")),
    pedantic = TRUE,
    quiet = TRUE
  )

3.0.2 Fit Model

Show the code
# number of MCMC chains
n_chains <- 4

# Run sampler
fit_quantifier <- model_quantifier$sample(
  data = stan_data, 
  seed = 2023,
  chains = n_chains,
  parallel_chains = n_chains,
  iter_warmup = 500,
  iter_sampling = 1000,
  refresh = 500,
  thin = 1,
  init = .1,
  adapt_delta = .8
)
  
# save fit
fit_quantifier$save_object(file =  here("fits", paste0(model_name_quantifier, "_custom_fit.RDS")))
Show the code
# load  fit
fit_quantifier <- readRDS(file =  here("fits", paste0(model_name_quantifier, "_custom_fit.RDS")))
Show the code
parameters <- parameters[parameters != "b_loc"]

3.0.3 Summary

Show the code
estimates_summary_quantifier <- fit_quantifier$summary(parameters)

3.0.4 Sampler Diagnostics

Show the code
# sampler diagnostics
fit_quantifier$sampler_diagnostics(format = "df") %>% 
  psych::describe(quant = c(.05,.95),) %>%
  round(2) %>%  
  as.data.frame() %>% 
  dplyr::select(median, min, Q0.05, Q0.95,  max) %>% 
  .[-c(7:9),]
              median    min  Q0.05  Q0.95    max
treedepth__     6.00   6.00   6.00   7.00   7.00
divergent__     0.00   0.00   0.00   0.00   0.00
energy__      629.85 500.88 567.93 696.98 766.12
accept_stat__   0.95   0.32   0.75   1.00   1.00
stepsize__      0.06   0.05   0.05   0.06   0.06
n_leapfrog__   63.00  63.00  63.00 127.00 255.00
Show the code
# convergence diagnostics
convergence_summary <- 
  fit_quantifier$draws(format = "df", variables = parameters) %>%
  summarise_draws(.x = ., "rhat", "ess_bulk", "ess_tail") %>%
  remove_missing() %>%
  select(-variable) %>%
  psych::describe(., quant = c(.05, .95)) %>%
  as.data.frame() %>%
  select(median, Q0.05, Q0.95, min, max)

convergence_summary %>% round(3)
           median    Q0.05    Q0.95     min      max
rhat        1.001    1.000    1.004   0.999    1.010
ess_bulk 1485.254  501.633 5409.092 297.833 7534.175
ess_tail 2355.948 1130.961 3162.272 619.152 3466.958

3.0.5 Effective sample size (ESS) & Rhat Plots

Show the code
# color scheme
color_scheme_set(scheme = "purple")

# Effective sample sizes
plot_neff <-
  mcmc_neff_hist(bayesplot::neff_ratio(fit_quantifier, pars = parameters),
                 binwidth = .01) +
  labs(title = "A") +
  guides(color = FALSE, fill = FALSE) +
  theme(
    legend.text = element_blank(),
    legend.key = element_blank(),
    title = element_text(size = 16, face = "bold")
  )
# Rhat
plot_rhat <-
  bayesplot::mcmc_rhat_hist(bayesplot::rhat(fit_quantifier, pars = parameters)) +
  labs(title = "B") +
  guides(color = FALSE, fill = FALSE) +
  theme(
    legend.text = element_blank(),
    legend.key = element_blank(),
    title = element_text(size = 16, face = "bold")
  ) +
  yaxis_text(on = TRUE)

# Combined plot
plot_diagnostics <- gridExtra::grid.arrange(plot_neff, plot_rhat, ncol = 2)

3.1 Parameter Plots

Show the code
# color scheme
color_scheme_set(scheme = "purple")

3.1.1 Person Competence: Location

Show the code
# order estimates by size
E_loc_ordered <-
  fit_quantifier$summary("E_loc") %>% 
  as.data.frame() %>% 
  arrange(median) %>% 
  select(variable) %>% 
  unlist()

# plot
  mcmc_intervals(fit_quantifier$draws(E_loc_ordered), point_est = "median", transformations = "log") +
  labs(
    subtitle = "Person Competence: Location",
    x = expression(log(E[i])),
    y = "Respondent"
  ) +
  scale_y_discrete(labels = NULL, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

3.1.2 Person Competence: Width

Show the code
# order estimates by size of theta
E_wid_ordered <- 
  str_replace(E_loc_ordered, "E_loc", "E_wid")

# plot

  mcmc_intervals(fit_quantifier$draws(E_wid_ordered), point_est = "median",transformations = "log") +
  labs(subtitle = "Person Competence: Width",
       x = expression(log(E_wid[i])),
       y = "Respondent") +
  scale_y_discrete(labels = NULL, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

3.1.3 Correlation: Person Competence Location vs. Width

Show the code
mcmc_intervals(fit_quantifier$draws("rho_E"), point_est = "median") +
  labs(
    title = "B",
    subtitle = "Correlation: Person Competence Location vs. Width",
    x = expression(Omega[i]),
    y = "Respondent"
  ) +
  xlim(-1,1)  +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

Show the code
fit_quantifier$draws("rho_E") %>% 
  bayestestR::describe_posterior(ci_method = "hdi")
Summary of Posterior Distribution

Parameter | Median |       95% CI |   pd |          ROPE | % in ROPE
--------------------------------------------------------------------
rho_E     |   0.63 | [0.51, 0.73] | 100% | [-0.10, 0.10] |        0%

3.1.4 Person Scaling Bias: Location

Show the code
# order estimates by size
a_loc_ordered <-
  fit_quantifier$summary("a_loc") %>% 
  as.data.frame() %>% arrange(median) %>% 
  select(variable) %>% 
  unlist()

# plot
mcmc_intervals(fit_quantifier$draws(a_loc_ordered), point_est = "median", transformations = "log") +
  labs(
    subtitle = "Person Scaling Bias: Location",
    x = expression(log(a[i])),
    y = "Respondent"
  ) +
  scale_y_discrete(labels = NULL, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

3.1.5 Person Shifting Bias: Width

Show the code
mcmc_intervals(fit_quantifier$draws("b_wid"), point_est = "median") +
  labs(
    subtitle = "Person Shifting Bias: Width",
    x = expression(b_wid[i]),
    y = "Respondent"
  ) +
  scale_y_discrete(labels = NULL, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

3.1.6 Item Discernability: Location

Show the code
 mcmc_intervals(fit_quantifier$draws("lambda_loc"), point_est = "median", transformations = "log") +
  labs(
    subtitle = "Item Discernability: Location",
    x = expression(log(lambda[loc])),
    y = "Item"
  ) +
  scale_y_discrete(labels = item_names_quantifier, expand = expansion(.02))  +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

3.1.7 Item Discernability: Width

Show the code
mcmc_intervals(fit_quantifier$draws("lambda_wid"), point_est = "median",transformations = "log") +
  labs(
    subtitle = "Item Discernability: Width",
    x = expression(log(lambda[wid])),
    y = "Item"
  ) +
  scale_y_discrete(labels = item_names_quantifier, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

3.1.8 Correlation: Item Discernibility Location vs. Width

Show the code
mcmc_intervals(fit_quantifier$draws("rho_lambda"), point_est = "median") +
  labs(
    x = expression(Omega[j])
  ) +
  xlim(-1,1) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

Show the code
fit_quantifier$draws("rho_lambda") %>% 
  bayestestR::describe_posterior(ci_method = "hdi")
Summary of Posterior Distribution

Parameter  | Median |         95% CI |     pd |          ROPE | % in ROPE
-------------------------------------------------------------------------
rho_lambda |  -0.47 | [-0.77, -0.06] | 98.00% | [-0.10, 0.10] |     2.68%

3.1.9 Correlation Between Location and Width

Show the code
  mcmc_intervals(fit_quantifier$draws(c("omega")), point_est = "median") +
  labs(
    subtitle = "Residual Correlation",
    x = expression(omega[j]),
    y = "Item"
  ) +
  scale_x_continuous(limits = c(-1,1)) +
  scale_y_discrete(labels = item_names_quantifier, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

3.1.10 Latent Truth Intervals

Show the code
consensus <- data.frame(
  idx = 1:J,
  name = item_names_quantifier,
  lower_consensus = fit_quantifier$summary("Tr_splx") %>%
    as.data.frame() %>%
    pull(median) %>%
    .[1:J],
  wid_consensus = fit_quantifier$summary("Tr_splx") %>%
    as.data.frame() %>%
    pull(median) %>%
    .[(J + 1):(J * 2)]
) %>%
  mutate(
    loc_consensus = lower_consensus + wid_consensus / 2,
    upper_consensus = lower_consensus + wid_consensus
  ) %>%
  arrange(lower_consensus)

consensus %>%
  # plot
  ggplot(aes(x = loc_consensus, y = 1:J)) +
  geom_errorbarh(aes(xmin = lower_consensus, xmax = upper_consensus), height = 0.3) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = seq(0, 1, .25),
    expand = expansion()
  ) +
  scale_y_continuous(breaks = 1:J, labels = item_names_quantifier[consensus$idx]) +
  labs(x = "Consensus Probabilities", y = NULL) +
  theme_icm(base_size = 14) +
  theme(
    panel.grid.x = element_blank(),
    
    axis.line.y = element_blank(),
    axis.ticks.x = element_line(colour = "#6d6d6e")
  )

Show the code
# save plot
ggsave(
  here("plots", "empirical_example", "consensus_intervals_custom_model.pdf"),
  width = 15,
  height = 8,
  units = "cm",
  scale = 1.3
)
Show the code
consensus %>% 
  mutate(across(c(lower_consensus, wid_consensus, loc_consensus, upper_consensus), ~round(., 3)*100)) %>% 
  select(-idx,-wid_consensus)
                 name lower_consensus loc_consensus upper_consensus
1               never             0.6           1.4             2.2
2        almost never             2.7           7.4            12.1
3              hardly             4.7          11.5            18.2
4              rarely             5.0          12.6            20.1
5        now and then            18.4          29.0            39.5
6        occasionally            18.6          29.6            40.5
7            possibly            19.7          33.9            48.1
8           sometimes            22.7          34.7            46.6
9         potentially            28.6          46.2            63.9
10 fifty-fifty chance            48.4          50.1            51.7
11         frequently            63.7          75.9            88.0
12              often            66.5          77.6            88.7
13   most of the time            70.3          81.2            92.1
14         constantly            78.0          87.2            96.4
15      almost always            84.4          90.8            97.2
16             always            96.9          98.2            99.4

4 Density Plots

4.0.1 Prepare Data

Show the code
df_long <- df_long %>%
  group_by(jj) %>%
  mutate(
    loc_mean = mean((x_L + x_U) / 2, na.rm = TRUE),
    wid_mean = mean(x_U - x_L, na.rm = TRUE),
    lower_mean = loc_mean - wid_mean / 2,
    upper_mean = loc_mean + wid_mean / 2,
    loc_mean_logit = mean(x_bvn_1, na.rm = TRUE),
    wid_mean_logit = mean(x_bvn_2, na.rm = TRUE),
    loc_median_logit = median(x_bvn_1, na.rm = TRUE),
    wid_median_logit = median(x_bvn_2, na.rm = TRUE)
  ) %>%
  ungroup()

splx <- bvn_to_simplex(df_long[, c("loc_mean_logit", "wid_mean_logit")])
df_long$lower_mean_logit <- splx[, 1]
df_long$upper_mean_logit <- 1 - splx[, 3]

splx <- bvn_to_simplex(df_long[, c("loc_median_logit", "wid_median_logit")])
df_long$lower_median_logit <- splx[, 1]
df_long$upper_median_logit <- 1 - splx[, 3]

rm(splx)

4.0.2 Plot: Verbal Quantifiers All

Show the code
df <- df_long %>%
  select(jj,
         name_en,
         x_L,
         x_U,
         lower_mean,
         upper_mean,
         lower_mean_logit,
         upper_mean_logit,
         lower_median_logit,
         upper_median_logit,
         truth) %>%
  remove_missing(vars = c("x_L", "x_U")) %>%
  full_join(., consensus, by = c("jj" = "idx"))

plot_density <-
  plot_intvls_aggregated(
    lower = df$x_L,
    upper = df$x_U,
    item_id = as.double(df$jj),
    lower_mean = df$lower_consensus * 100,
    upper_mean = df$upper_consensus * 100,
    lower_mean_logit = df$lower_median_logit * 100,
    upper_mean_logit = df$upper_median_logit * 100,
    item_name = df$name_en,
    facet_wrap = TRUE,
    min = 0,
    max = 100,
    step_size = 1,
    show_quantiles = FALSE
  )

plot_density

Show the code
# save plot
ggsave(
  plot = plot_density,
  here("plots", "empirical_example", "verbal_quantifier_densities_all.pdf"),
  width = 15,
  height = 15,
  units = "cm",
  scale = 1.7
)

4.0.3 Plot: Verbal Quantifiers Selection for Article

Show the code
df <- df_long %>%
  select(
    jj,
    name_en,
    x_L,
    x_U,
    lower_mean,
    upper_mean,
    lower_mean_logit,
    upper_mean_logit,
    lower_median_logit,
    upper_median_logit,
    truth
  ) %>%
  dplyr::filter(name_en %in% c("fifty-fifty chance", "hardly", "potentially", "never", "always")) %>%
  full_join(., consensus, by = c("jj" = "idx")) %>%
  remove_missing(vars = c("x_L", "x_U"))


# Layout for facet wrap
design <- matrix(c(2, 2,  4, 1, 3, 5), 3, 2, byrow = TRUE)

# plot
plot_selection <-
  plot_intvls_aggregated(
    lower = df$x_L,
    upper = df$x_U,
    item_id = as.double(df$jj),
    lower_mean = df$lower_consensus * 100,
    upper_mean = df$upper_consensus * 100,
    lower_mean_logit = df$lower_median_logit * 100,
    upper_mean_logit = df$upper_median_logit * 100,
    item_name = (df$name_en),
    #truth = df$truth,
    facet_wrap = TRUE,
    design = design,
    ncol = 2,
    min = 0,
    max = 100,
    step_size = 1,
    show_quantiles = FALSE
  )
plot_selection

Show the code
# save plot
ggsave(
  plot = plot_selection,
  here("plots", "empirical_example", "verbal_quantifier_densities_example.pdf"),
  width = 15,
  height = 10,
  units = "cm",
  scale = 1.6
)

4.0.4 Plot for slides

Show the code
df <- df_long %>%
    select(
    jj,
    name_en,
    x_L,
    x_U,
    lower_mean,
    upper_mean,
    lower_mean_logit,
    upper_mean_logit,
    truth
  ) %>%
  dplyr::filter(name_en %in% c("fifty-fifty chance")) %>%
  full_join(., consensus, by = c("jj" = "idx")) %>%
  remove_missing(vars = c("x_L", "x_U")) 


# plot
plot_selection <-  plot_intvls_aggregated(
  lower = df$x_L,
  upper = df$x_U,
  item_id = as.double(df$jj),
  lower_mean = df$lower_consensus*100,
  upper_mean = df$upper_consensus*100,
  lower_mean_logit = df$lower_mean_logit*100,
  upper_mean_logit = df$upper_mean_logit*100,
  item_name = (df$name_en),
  truth = df$truth,
  facet_wrap = FALSE,
  ncol = 2,
  min = 0,
  max = 100,
  step_size = 1,
  show_quantiles = FALSE
)
#plot_selection

# save plot
ggsave(
  plot = plot_selection,
  here("plots", "empirical_example", "verbal_quantifier_densities_50-50.pdf"),
  width = 15,
  height = 5,
  units = "cm",
  scale = 1.6
)

ggsave(
  plot = plot_selection,
  here("plots", "empirical_example", "verbal_quantifier_densities_50-50.svg"),
  width = 15,
  height = 5,
  units = "cm",
  scale = 1.6
)

4.0.5 Plot for slides

Show the code
df <- df_long %>%
  select(
    jj,
    name,
    name_en,
    x_L,
    x_U,
    lower_mean,
    upper_mean,
    lower_mean_logit,
    upper_mean_logit,
    truth
  ) %>%
  dplyr::filter(name %in% c("meistens")) %>%
  full_join(., consensus, by = c("jj" = "idx")) %>%
  remove_missing(vars = c("x_L", "x_U")) 


# plot
plot_selection <-  plot_intvls_aggregated(
  lower = df$x_L,
  upper = df$x_U,
  item_id = as.double(df$jj),
  lower_mean = df$lower_consensus*100,
  upper_mean = df$upper_consensus*100,
  lower_mean_logit = df$lower_mean_logit*100,
  upper_mean_logit = df$upper_mean_logit*100,
  item_name = (df$name_en),
  truth = df$truth,
  facet_wrap = FALSE,
  ncol = 2,
  min = 0,
  max = 100,
  step_size = 1,
  show_quantiles = FALSE
)
#plot_selection

# save plot
ggsave(
  plot = plot_selection,
  here("plots", "empirical_example", "verbal_quantifier_densities_mostly.pdf"),
  width = 15,
  height = 5,
  units = "cm",
  scale = 1.6
)

ggsave(
  plot = plot_selection,
  here("plots", "empirical_example", "verbal_quantifier_densities_mostly.svg"),
  width = 15,
  height = 5,
  units = "cm",
  scale = 1.6
)

4.1 Responses vs. Proficiency

Show the code
# prepare data
proficiency_loc <- 
  fit_quantifier$summary("E_loc") %>% as.data.frame() %>% select(median) %>% 
  rename(proficiency_loc = median) %>%
  mutate(ii = 1:I,
         proficiency_loc = scale(log(proficiency_loc)))

proficiency_wid <- 
  fit_quantifier$summary("E_wid") %>% as.data.frame() %>% select(median) %>%
  rename(proficiency_wid = median) %>%
  mutate(ii = 1:I,
         proficiency_wid = scale(log(proficiency_wid))) 

proficiency_medians <-
  full_join(proficiency_loc, proficiency_wid, by = c("ii" = "ii")) %>%
  full_join(., df_long %>%
              select(ii, x_L_u, x_U_u, name), by = c("ii" = "ii")) %>%
  mutate(
    proficiency = (proficiency_loc + proficiency_wid) / 2,
    ii_ranked = factor(proficiency) %>% as.numeric()
  ) %>%
  arrange(proficiency) %>%
  dplyr::filter(name == "Fuenfzig-Fuenfzig Chance")


# plot
plot_proficiency_vs_intervals <- 
  proficiency_medians %>%
  ggplot() +
  geom_rect(aes(
    xmin = 40,
    xmax = 60,
    ymin = 0,
    ymax = max(ii_ranked) + 1
  ), color = "grey80") +
  geom_errorbar(aes(
    y = (ii_ranked),
    xmin = x_L_u,
    xmax = x_U_u
  ),  
  alpha = .5,
  width = 2,
  linewidth = .5) +
  geom_point(
    aes(x = pnorm(proficiency), y = ii_ranked), 
    shape = 16,
    col = ggokabeito::palette_okabe_ito()[5]
    ) +
  scale_x_continuous(
    limits = c(0, 1),
    breaks = seq(0, 1, .25),
    labels = c("0", ".25", ".50", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(expand = expansion(0, .2)) +
  labs(
    x = "Response Interval / <span style='color: #0072B2;'>Transformed Proficiency</span>", 
    y = "Respondent") +
  coord_cartesian(clip = "off") +
  theme_icm() +
  theme(
    axis.title.x = ggtext::element_markdown(),
    axis.line = element_line(colour = "#6d6d6e", size = .3),
    axis.ticks.x = element_line(colour = "#6d6d6e", size = .3),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    panel.grid = ggplot2::element_blank(
    )
  )

ggsave(
  plot = plot_proficiency_vs_intervals,
  filename = here("plots", "empirical_example", "proficiency_vs_intervals.pdf"),
  width = 15,
  height = 10,
  units = "cm",
  scale = 1.4
)

plot_proficiency_vs_intervals

5 Prior vs. Posterior

5.1 True Intervals

Show the code
# get posterior
consensus <- data.frame(
  Tr_loc_splx_post = fit_quantifier$draws("Tr_loc_splx") %>%
    unlist() %>%
    as.vector(),
  Tr_wid_splx_post = fit_quantifier$draws("Tr_wid_splx") %>%
    unlist() %>%
    as.vector()
)  %>%
  mutate(jj = (rep(1:J, each = nrow(.) / J))) %>%
  group_by(jj)  %>%
  mutate(
    text_x = mean(Tr_loc_splx_post),
    text_y_mean = mean(Tr_wid_splx_post),
    text_y_max = max(Tr_wid_splx_post)
  ) %>%
  ungroup()  %>%
  full_join(., df_long %>%
              select(jj, name_en), by = c("jj" = "jj")) |> 
  slice_sample(n = 1000, by = jj) 

# compute prior samples
n_samples <- 1e6
Tr_prior_samples <- data.frame(Tr_wid_prior = rbeta(n_samples, 1.2, 3))
Tr_prior_samples$Tr_loc_prior <-
  rbeta(n_samples, 1, 1) * (1 - Tr_prior_samples$Tr_wid_prior) + Tr_prior_samples$Tr_wid_prior / 2

# Plot
consensus %>%
  ggplot() +
  # prior density
  geom_density_2d_filled(
    data = Tr_prior_samples,
    aes(
      x = Tr_loc_prior,
      y = Tr_wid_prior,
      fill = after_stat(as.numeric(level) / max(as.numeric(level)))
    ),
    bins = 100,
    # Increase number of bins
    alpha = 0.5,
    show.legend = TRUE
  ) +
  # legend for prior density
  scale_fill_viridis_c(name = "Prior Density", alpha = .5, 
                       breaks = c(0,.25,.5,.75,1), 
                       labels = c("0",".25",".5",".75","1")) +
  # mask everything outside of the ternary space
  annotate(
    geom = "polygon",
    x = c(0,-Inf, -Inf,  .5),
    y = c(0,  0, Inf, Inf),
    color = "white",
    fill = "white"
  ) +
  annotate(
    geom = "polygon",
    x = c(1, Inf, Inf, .5),
    y = c(0, 0, Inf, Inf),
    color = "white",
    fill = "white"
  ) +
  # triangle lines
  geom_segment(aes(
    x = 0,
    y = 0,
    xend = .5,
    yend = 1
  )  ,
  col = "#6d6d6e",
  lineend = "round") +
  geom_segment(aes(
    x = 1,
    y = 0,
    xend = .5,
    yend = 1
  ),
  col = "#6d6d6e",
  lineend = "round") +
  # posterior samples
  geom_point(
    data = consensus,
    aes(x = Tr_loc_splx_post, y = Tr_wid_splx_post, col = factor(jj)),
    size = .3,
    alpha = .05,
    shape = 16
  ) +
  geom_segment(
    data = consensus |> distinct(jj, .keep_all = TRUE),
    aes(
    x = text_x,
    y = text_y_mean,
    xend = text_x,
    yend = text_y_max + .002
  )) +
  geom_label(
    data = consensus |> distinct(jj, .keep_all = TRUE),
    aes(x = text_x, y  = text_y_max, label = name_en),
    label.size = 0,
    alpha = .5,
    label.padding = unit(0.25, "lines"),
    nudge_y = .02,
    size.unit = "pt", 
    size = 12,
  ) +
  # format scales
  scale_x_continuous(
    limits = c(-.05, 1.05),
    breaks = seq(0, 1, .25),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  coord_equal(clip = "off") +
  labs(x = "Location", y = "Width") +
  # suppress legends
  guides(col = FALSE, fill = FALSE) +
  theme_icm() +
  theme(
    panel.grid = ggplot2::element_line(colour = "white"),
    plot.margin = margin(.3, 0, -.03, 0, "cm"),
    legend.position = "right",
    legend.title = element_text(
      margin = margin(b = 15),
      # Top, Right, Bottom, Left
      size = 14
      )
  )

5.1.1 Plot for article with selected quantifiers

Show the code
# Data
consensus_selection <- 
  consensus |> 
  dplyr:: filter(
    name_en %in% c(
      "never",
      "rarely",
      "now and then",
      "possibly",
      "potentially",
      "fifty-fifty chance",
      "frequently",
      "almost always",
      "always"
    )
  ) 
# Plot
ggplot() +
  # prior density
  geom_density_2d_filled(
    data = Tr_prior_samples,
    aes(
      x = Tr_loc_prior,
      y = Tr_wid_prior,
      fill = after_stat(as.numeric(level) / max(as.numeric(level)))
    ),
    bins = 100,
    # Increase number of bins
    alpha = 0.5,
    show.legend = TRUE
  ) +
  # legend for prior density
  scale_fill_viridis_c(name = "Prior Density", alpha = .5, 
                       breaks = c(0,.25,.5,.75,1), 
                       labels = c("0",".25",".5",".75","1")) +
  # mask everything outside of the ternary space
  annotate(
    geom = "polygon",
    x = c(0, -Inf, -Inf, .5),
    y = c(0, 0, Inf, Inf),
    color = "white",
    
    fill = "white"
  ) +
  annotate(geom = "polygon",
           x = c(1, Inf, Inf, .5),
           y = c(0, 0,  Inf, Inf),
    color = "white",
    fill = "white"
  ) +
  # triangle lines
  geom_segment(aes(
    x = 0,
    y = 0,
    xend = .5,
    yend = 1
  )  ,
  col = "#6d6d6e",
  lineend = "round") +
  geom_segment(aes(
    x = 1,
    y = 0,
    xend = .5,
    yend = 1
  ),
  col = "#6d6d6e",
  lineend = "round") +
  # posterior samples
  geom_point(
    data = consensus_selection,
    aes(x = Tr_loc_splx_post, y = Tr_wid_splx_post),
    col = "#d95f02",
    size = .3,
    alpha = .3,
    shape = 16
  ) +
  geom_segment(
    data = consensus_selection |> distinct(jj, .keep_all = TRUE),
    aes(
    x = text_x,
    y = text_y_mean,
    xend = text_x,
    yend = text_y_max + .002
  )) +
  geom_label(
    data = consensus_selection |> distinct(jj, .keep_all = TRUE),
    aes(x = text_x, y  = text_y_max, label = name_en),
    label.size = 0,
    alpha = .5,
    label.padding = unit(0.25, "lines"),
    nudge_y = .02,
    size.unit = "pt", 
    size = 12,
  ) +
  # format scales
  scale_x_continuous(
    limits = c(-.05, 1.05),
    breaks = seq(0, 1, .25),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  coord_equal() +
  labs(x = "Location", y = "Width") +
  # suppress legends
  guides(col = FALSE) +
  theme_icm() +
  theme(
    panel.grid = ggplot2::element_line(colour = "white"),
    plot.margin = margin(.3, 0, -.03, 0, "cm"),
    legend.position = "right",
    legend.title = element_text(
      margin = margin(b = 15),
      # Top, Right, Bottom, Left
      size = 14
      )
  )

Show the code
ggsave(
  plot = last_plot(),
  filename = here("plots", "empirical_example", "prior_vs_posterior_quantifiers_selection.pdf"),
  width = 15,
  height = 10,
  units = "cm",
  scale = 1.4
)

ggsave(
  plot = last_plot(),
  filename = here("plots", "empirical_example", "prior_vs_posterior_quantifiers_selection.svg"),
  width = 15,
  height = 10,
  units = "cm",
  scale = 1.4
)

5.2 Hyper-Priors

5.2.1 sigma_lambda: Location

Show the code
colors <- c("Prior" = "lightblue", "Posterior" = "purple")

sigma_lambda_loc <- data.frame(posterior = fit_quantifier$draws("sigma_lambda[1]", format = "list") %>% unlist() %>%  as.vector())
sigma_lambda_loc$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_lambda_loc %>%
  ggplot() +
  geom_density(aes(prior, fill = "Prior"), alpha = 1) +
  geom_density(aes(exp(posterior * .5 + log(.5)), fill = "Posterior"), alpha = .3) +
  scale_x_continuous(limits = c(NA, 3), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  scale_fill_manual(values = colors) +
  labs(x = "Prior / Posterior", y = "Density", fill = NULL)  +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank(), 
        legend.position = "right")

5.2.2 sigma_lambda: Width

Show the code
sigma_lambda_wid <- data.frame(
  posterior = fit_quantifier$draws("sigma_lambda[2]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_lambda_wid$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_lambda_wid %>% 
  ggplot() +
  geom_density(
    aes(prior,fill = "Prior"),
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5)), fill = "Posterior"),
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  scale_fill_manual(values = colors) +
  labs(x = "Prior / Posterior", y = "Density", fill = NULL) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank(), 
        legend.position = "right")

5.2.3 sigma_E: Location

Show the code
sigma_E_loc <- data.frame(
  posterior = fit_quantifier$draws("sigma_I[1]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_E_loc$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_E_loc %>% 
  ggplot() +
  geom_density(
    aes(prior,fill = "Prior"),
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5)), fill = "Posterior"),
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  scale_fill_manual(values = colors) +
  labs(x = "Prior / Posterior", y = "Density", fill = NULL) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank(), 
        legend.position = "right")

5.2.4 sigma_E: Width

Show the code
sigma_E_wid <- data.frame(
  posterior = fit_quantifier$draws("sigma_I[2]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_E_wid$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_E_wid %>% 
  ggplot() +
  geom_density(
    aes(prior, fill = "Prior"),
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5)), fill = "Posterior"),
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  scale_fill_manual(values = colors) +
  labs(x = "Prior / Posterior", y = "Density", fill = NULL) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank(), 
        legend.position = "right")

5.2.5 mu_E: Location

Show the code
mu_E_loc <- data.frame(
  posterior = fit_quantifier$draws("mu_E[1]", format = "list") %>% unlist() %>%  as.vector()
)
mu_E_loc$prior <- rnorm(nrow(sigma_lambda_loc))

mu_E_loc %>% 
  ggplot() +
  geom_density(
    aes(prior, fill = "Prior"),
    alpha = 1  ) +
  geom_density(aes(posterior, fill = "Posterior"),
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(-4, 4),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  scale_fill_manual(values = colors) +
  labs(x = "Prior / Posterior", y = "Density", fill = NULL) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank(), 
        legend.position = "right")

5.2.6 mu_E: Width

Show the code
mu_E_wid <- data.frame(
  posterior = fit_quantifier$draws("mu_E[2]", format = "list") %>% unlist() %>%  as.vector()
)
mu_E_wid$prior <-  rnorm(nrow(sigma_lambda_loc))

mu_E_wid %>% 
  ggplot() +
  geom_density(
    aes(prior, fill = "Prior"),
    alpha = 1  ) +
  geom_density(aes(posterior, fill = "Posterior"),
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(-4, 4),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  scale_fill_manual(values = colors) +
  labs(x = "Prior / Posterior", y = "Density", fill = NULL) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank(), 
        legend.position = "right")

6 Posterior Predictive Checks

Show the code
sample <- sample(1:1000, 100)
Y_ppc_loc <- fit_quantifier$draws("Y_ppc_loc", format = "matrix")[sample,]
Y_ppc_wid <- fit_quantifier$draws("Y_ppc_wid", format = "matrix")[sample,]
Y_ppc_loc_splx <- fit_quantifier$draws("Y_ppc_loc_splx", format = "matrix")[sample,]
Y_ppc_wid_splx <- fit_quantifier$draws("Y_ppc_wid_splx", format = "matrix")[sample,]

6.0.1 Location: Bivariate Normal

Show the code
bayesplot::ppc_dens_overlay(y = df_long$x_bvn_1, yrep = Y_ppc_loc) +
  xlim(-10, 10)

6.0.2 Width: Bivariate Normal

Show the code
bayesplot::ppc_dens_overlay(
  y = df_long$x_bvn_2,
  yrep = Y_ppc_wid
  ) +
  xlim(-10,10) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

6.0.3 Location: Bounded

Show the code
bayesplot::ppc_dens_overlay(y = df_long$x_M_u, yrep = Y_ppc_loc_splx) +
  xlim(0, 1)

6.0.4 Width: Bounded

Show the code
bayesplot::ppc_dens_overlay(y = df_long$x_W_u, yrep = Y_ppc_wid_splx) +
  xlim(0, 1) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

6.0.5 Lower Interval Bound

Show the code
bayesplot::ppc_dens_overlay(y = df_long$x_L_u, yrep = Y_ppc_loc_splx - Y_ppc_wid_splx/2) +
  xlim(0, 1)

6.0.6 Upper Interval Bound

Show the code
bayesplot::ppc_dens_overlay(y = df_long$x_U_u, yrep = Y_ppc_loc_splx + Y_ppc_wid_splx/2) +
  xlim(0, 1) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

7 Re-Fit Without Control Items

7.1 Stan Data

Show the code
### Stan data declaration
df_long_no_controls <- 
  df_long %>% 
  dplyr::filter(!name %in% c("Fuenfzig-Fuenfzig Chance", "immer", "niemals")) %>% 
  mutate(
    ii = as.integer(as.factor(ii)),
    jj = as.integer(as.factor(jj))
  )

## stan data list
stan_data_no_controls <- list(
  I = length(unique(df_long_no_controls$ii)),
  J = length(unique(df_long_no_controls$jj)),
  N = nrow(df_long_no_controls),
  ii =  df_long_no_controls$ii,
  jj = df_long_no_controls$jj,
  nn = c(1:nrow(df_long_no_controls)),
  padding = .01,
  Y_splx = cbind(
    df_long_no_controls$x_splx_1,
    df_long_no_controls$x_splx_2,
    df_long_no_controls$x_splx_3
  ) %>%
    as.matrix()
)

7.1.1 Fit Model

Show the code
# number of MCMC chains
n_chains <- 4

# Run sampler
fit_quantifier_no_controls <- model_quantifier$sample(
  data = stan_data_no_controls, 
  seed = 2023,
  chains = n_chains,
  parallel_chains = n_chains,
  iter_warmup = 500,
  iter_sampling = 1000,
  refresh = 500,
  thin = 1,
  init = .1,
  adapt_delta = .8
)
  
# save fit
fit_quantifier_no_controls$save_object(file =  here("fits", paste0(model_name_quantifier, "_no_controls_fit.RDS")))
Show the code
# load  fit
fit_quantifier_no_controls <- readRDS(file =  here("fits", paste0(model_name_quantifier, "_no_controls_fit.RDS")))

7.1.2 Item Discernability: Location

Show the code
mcmc_intervals(
  fit_quantifier_no_controls$draws("lambda_loc"),
  point_est = "median",
  transformations = "log"
) +
  labs(subtitle = "Item Discernability: Location",
       x = expression(log(lambda[loc])),
       y = "Item") +
  scale_y_discrete(labels = item_names_quantifier, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

7.1.3 Item Discernability: Width

Show the code
mcmc_intervals(
  fit_quantifier_no_controls$draws("lambda_wid"),
  point_est = "median",
  transformations = "log"
) +
  labs(subtitle = "Item Discernability: Width",
       x = expression(log(lambda[wid])),
       y = "Item") +
  scale_y_discrete(labels = item_names_quantifier, expand = expansion(.02)) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

7.2 Discernibility: Correlation Between Location and Width

Show the code
mcmc_intervals(fit_quantifier_no_controls$draws("rho_lambda"), point_est = "median") +
  labs(
    x = expression(Omega[j])
  ) +
  xlim(-1,1) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

Show the code
fit_quantifier_no_controls$draws("rho_lambda") %>% 
  bayestestR::describe_posterior(ci_method = "hdi")
Summary of Posterior Distribution

Parameter  | Median |       95% CI |     pd |          ROPE | % in ROPE
-----------------------------------------------------------------------
rho_lambda |   0.80 | [0.43, 0.98] | 99.90% | [-0.10, 0.10] |        0%

7.3 Proficiency: Correlation Between Location and Width

Show the code
mcmc_intervals(fit_quantifier_no_controls$draws("rho_E"), point_est = "median") +
  labs(
    x = expression(Omega[j])
  ) +
  xlim(-1,1) +
  theme_icm(base_size = 16) +
  theme(panel.grid = element_blank())

Show the code
fit_quantifier_no_controls$draws("rho_E") %>% 
  bayestestR::describe_posterior(ci_method = "hdi")
Summary of Posterior Distribution

Parameter | Median |       95% CI |   pd |          ROPE | % in ROPE
--------------------------------------------------------------------
rho_E     |   0.50 | [0.34, 0.64] | 100% | [-0.10, 0.10] |        0%

8 Session Info

Show the code
sessionInfo()
R version 4.5.0 (2025-04-11 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=German_Germany.utf8  LC_CTYPE=German_Germany.utf8   
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.utf8    

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] showtext_0.9-7    showtextdb_3.0    sysfonts_0.8.9    extraDistr_1.10.0
 [5] ggtext_0.1.2      ggh4x_0.3.1       ggokabeito_0.1.0  ggpubr_0.6.0     
 [9] osfr_0.2.9        here_1.0.1        psych_2.5.6       svglite_2.2.1    
[13] cowplot_1.1.3     rmarkdown_2.29    bayestestR_0.16.0 posterior_1.6.1  
[17] bayesplot_1.13.0  cmdstanr_0.8.0    readxl_1.4.5      yardstick_1.3.2  
[21] lubridate_1.9.4   forcats_1.0.0     stringr_1.5.1     dplyr_1.1.4      
[25] purrr_1.0.4       readr_2.1.5       tidyr_1.3.1       tibble_3.3.0     
[29] ggplot2_3.5.2     tidyverse_2.0.0   quarto_1.4.4      pacman_0.5.1     

loaded via a namespace (and not attached):
 [1] mnormt_2.1.1         gridExtra_2.3        rlang_1.1.6         
 [4] magrittr_2.0.3       matrixStats_1.5.0    compiler_4.5.0      
 [7] reshape2_1.4.4       systemfonts_1.2.3    vctrs_0.6.5         
[10] httpcode_0.3.0       pkgconfig_2.0.3      fastmap_1.2.0       
[13] backports_1.5.0      labeling_0.4.3       markdown_2.0        
[16] tzdb_0.5.0           ps_1.9.1             ragg_1.4.0          
[19] xfun_0.52            cachem_1.1.0         litedown_0.7        
[22] jsonlite_2.0.0       later_1.4.2          broom_1.0.8         
[25] parallel_4.5.0       R6_2.6.1             stringi_1.8.7       
[28] RColorBrewer_1.1-3   car_3.1-3            cellranger_1.1.0    
[31] Rcpp_1.0.14          knitr_1.50           timechange_0.3.0    
[34] tidyselect_1.2.1     rstudioapi_0.17.1    abind_1.4-8         
[37] yaml_2.3.10          curl_6.4.0           processx_3.8.6      
[40] plyr_1.8.9           lattice_0.22-7       withr_3.0.2         
[43] evaluate_1.0.4       isoband_0.2.7        xml2_1.3.8          
[46] pillar_1.10.2        carData_3.0-5        tensorA_0.36.2.1    
[49] checkmate_2.3.2      insight_1.3.0        distributional_0.5.0
[52] generics_0.1.4       rprojroot_2.0.4      hms_1.1.3           
[55] commonmark_1.9.5     scales_1.4.0         glue_1.8.0          
[58] tools_4.5.0          ggsignif_0.6.4       fs_1.6.6            
[61] grid_4.5.0           datawizard_1.1.0     nlme_3.1-168        
[64] Formula_1.2-5        cli_3.6.5            textshaping_1.0.1   
[67] viridisLite_0.4.2    gtable_0.3.6         rstatix_0.7.2       
[70] digest_0.6.37        crul_1.5.0           htmlwidgets_1.6.4   
[73] farver_2.1.2         memoise_2.0.1        htmltools_0.5.8.1   
[76] lifecycle_1.0.4      httr_1.4.7           MASS_7.3-65         
[79] gridtext_0.1.5